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Abstract 

We consider a minisuperspace model for a closed universe with small and 
positive cosmological constant A, filled with a massive scalar field conformally 
coupled to gravity. In the quantum version of this model, the universe may 
undergo a tunneling transition through an effective barrier between regions 
of small and large scale factor. We solve numerically the minisuperspace 
Wheeler-De Witt equation with tunneling boundary conditions for the wave 
function of the universe, and find that tunneling in quantum cosmology is 
quite different from that in quantum mechanics. Namely, the matter degree 
of freedom gets excited under the barrier, provided its interaction with the 
scale factor is not too weak, and makes a strong back reaction onto tunneling. 
In the semiclassical limit of small A, the matter energy behind the barrier is 
close to the height of the barrier: the system "climbs up" the barrier, and then 
evolves classically from its top. These features are even more pronounced for 
inhomogeneous modes of matter field. Extrapolating to field theory we thus 
argue that high momentum particles are copiously created in the tunneling 
process. Nevertheless, we find empirical evidence for the semiclassical-type 
scaling with A of the wave function under and behind the barrier. 
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1 Introduction and summary. 



A traditional formalism in quantum cosmology is based on the Wlieeler-De Witt 
equation [H Ej, often truncated to minisuperspace. Though this equation is not 
applicable at length scales near and below the Planck scale, it should correctly 
describe quantum gravity phenomena at larger distances. Among these phenomena, 
penetration through classically forbidden regions is of particular interest. Indeed, 
various proposals for "quantum creation of the Universe" jSl IH El 13 El invoke this 
possible quantum gravity effect in one or another way (for discussion and references 
see, e.g., Refs. PHUl). 

A prototype example is an empty closed Universe with small and positive cosmo- 
logical constant A. In the minisuperspace approximation, the only dynamical degree 
of freedom is the radius of the Universe a. In Planck units and with appropriate 
rescaling, the gravitational action has the following form (see, e.g., Refs. IHllin] and 
references therein). 



where t] is the time variable (in = 1 gauge, t] has the meaning of conformal time), 
TTa is momentum conjugate to a, A^ is the lapse function and 



is the (conformal) Hamiltonian for the scale factor. Correspondingly, the minisu- 
perspace Wheeler-De Witt equation is 



where ^ is the wave function of the Universe, and a and Tia are treated as operators 
with standard commutation relation^ In the coordinate representation one has 
^ = \t'(a) and tTq = —id/ da. One may modify eq. (jT} slightly by adding a spatially 
homogeneous mode of a massless conformal scalar field. Then the Wheeler-De Witt 
equation becomes 



where e is the conformal energy of the scalar mode. 

There are various proposals for the wave function of the Universe |H1 ISl El (HI 
which, in the language of the Wheeler-De Witt equation, correspond to different 
boundary conditions for \l/(a). In this paper we consider tunneling boundary con- 
ditions [3171 m • These resemble closely the boundary conditions corresponding to 
tunneling transitions in conventional quantum mechanics. Namely, eq. ^ has the 
form of the stationary Schrodinger equation at energy e in a potential 





= , 




+ = 



(2) 




The operator ordering ambiguity, inherent in eq. (Q, is unimportant for our purposes. 
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Figure 1: Potential barrier in eq. (0). 

see Fig. ^ At small A (in Planck units) the system is semiclassical: it is straightfor- 
ward to see that A plays the role of h. At sufficiently large e (again in Planck units), 
but e < Vmax = 1/(16A), there are classically allowed regions on the two sides of 
the barrier. In conventional quantum mechanics, tunneling from small a to large 
a would be described by the wave function which contains only outgoing waves as 
a +00; if there are degrees of freedom other than a, then the incoming component 
of the wave function has to be specified on the left of the barrier. In this paper we 
consider the same boundary conditions in the context of quantum cosmology. 

At small A, the decaying part of the wave function in the classically forbidden 
region is 

^ - / ^«V2(V(a) - e) 
W(a) oc e , (3) 

so the wave function behind the barrier is suppressed by the standard semiclassical 
exponential 

-J da ^j2{V{a) - e) 
*la^+oo~e -1 (4) 

(assuming that \1/ ~ 1 in the classically allowed region on the left of the barrier), 
where Oi and 02 are the two turning points. 

The analogy to the stationary Schrodinger equation is no longer perfect if matter 
degrees of freedom (or gravitons) are added. Consider, as an example, a massive 
scalar field conformally coupled to gravity, whose action is 
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In the minisuperspace approximation, the scalar field is spatially homogeneous, 

a{rj) ' 

and its action, with appropriate rescaling, is 

S<), = J drj {(pn^ - NH^) , 

where 

7^<^ = i + i(l + mV)0^ (5) 

We note in passing that if one blindly considered an inhomogeneous mode of con- 
formal momentum fc, then the scalar Hamiltonian would be 

= Y + ^(^' + i + "^W- (6) 

With the scalar field added, the Wheeler-De Witt equation becomes 

{na + n^)^ + = Q . (7) 

Unlike in conventional quantum mechanics, this equation contains a negative "ki- 
netic term" for the scale factor, (— 7r^/2), and the potential V{a) enters with negative 
sign also. This feature, which at first sight appears technical (and is often ignored in 
literature, see, however, Ref. J2]); in fact refiects the difference between tunneling in 
quantum cosmology and in quantum mechanics. Consider, as an example, tunneling 
from the ground state of the 0-oscillator. In quantum mechanics, the excitation of a 
fluctuating degree of freedom (similar to the fleld 0) would lower the energy of the 
tunneling degree of freedom (the scale factor a), and thus suppress tunneling even 
further. Hence, the fluctuating degree of freedom does not get excited much and 
the semiclassical exponent (jH) does not change (the fluctuating degree of freedom 
affects the pre-exponential factor only). In quantum cosmology, on the contrary, the 
excitation of the fleld makes the tunneling of the Universe exponentially easier jjj: 
crudely speaking, the effective value of e gets increased. Hence, the 0-quanta tend 
to be copiously created. Because of their back reaction, the semiclassical expansion 
about the solution Q may break down, and the entire picture of tunneling in quan- 
tum cosmology may get substantially modifled when matter degrees of freedom or 
gravitons are added. 

It is worth pointing out, that at given e there exists only one real Euclidean 
solution in the theory with the action {Sa + S",^). This solution, the instanton (pe- 
riodic instanton for e 7^ 0), has = and coincides with the instanton in a theory 
truncated to a single variable a. The above arguments imply that this instanton 
may become completely irrelevant if matter degrees of freedom and / or gravitons are 
kept. 

In this paper we substantiate our qualitative observations by a numerical study of 
quantum systems of the Wheeler-De Witt type. Namely, we consider systems with 
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two degrees of freedom whose own Hamiltonians are of opposite signs. To obtain a 
clear picture, we modify eq. ((Tj) in such a way that the potential V{a) is well localized 
in a, and the interaction between (p and a is also localized in a. Then the tunneling 
boundary conditions are well defined even for a finite interval of a (which one has 
to work with in a numerical study), while the peculiarities of the Wheeler-De Witt 
equation remain there. To verify that the main results of our study are not sensitive 
to the details of the a-0 interaction, we modify the interaction in different ways, and 
consider three different models in subsequent sections. Our particular systems are 
described in section 2.1. 

If the interaction between and a is weak enough, the matter degree of free- 
dom may be treated as perturbation (provided it is in its ground state on the left 
of the barrier, otherwise the parameter e should be redefined). To the zeroth ap- 
proximation, the wave function in the classically forbidden region is then given by 
eq. (jni). In section 2.2 we present a perturbative treatment of our systems about the 
solution (jni) and get a heuristic understanding of the range of parameters in which 
the perturbation theory breaks down. Unlike in conventional quantum mechanics, 
this range includes part of the would-be semiclassical region of small A. We will be 
mostly interested in this range of parameters. 

Our main purpose is to study our systems at the quantum level. We have solved 
the Wheeler-De Witt equation for these systems numerically, without making use 
of any approximation (except for approximations inevitable in numerical studies, 
such as discretization, truncation of the matter Hamiltonian to finite number of 
levels, etc.; these approximations are safely under control). We describe our results 
in section 3. We begin with the numerical analysis of the classical evolution (section 
3.1) to see that, in a wide range of parameters, the potential barrier cannot be 
overcome classically. We then proceed in section 3.2 to the quantum treatment. 
We find that unless the interaction between a and is weak, the matter degree of 
freedom does indeed get strongly excited under the barrier, and the system behind 
the barrier contains a large number of 0-quanta. Furthermore, in the region of the 
parameter space which in conventional quantum mechanics would be adequately 
described semiclassically, with the wave function (jH]), the energy of the matter degree 
of freedom is close, behind the barrier, to the height of the barrier Vmax ~ V^- 
The system "climbs up" the potential barrier, and the classically allowed region 
effectively starts again just right of the point of maximum, a^ax- These features are 
even more pronounced for a larger a-independent term in the oscillator frequency, k'^ 
in eq. (jH)), again in contrast to quantum mechanics. Even if the interaction between 
matter and gravity is weak (i.e., the parameter m in eq. © is small), the 0-quanta 
are copiously created due to tunneling, provided that their frequency is high {k is 
large enough). Extrapolating to field theory in the tunneling Universe, we thus 
argue that high momentum modes get strongly excited in the tunneling process. 

In section 3.3 we concentrate on scaling properties in A of our quantum solutions 
in the region of small A. We find, somewhat surprisingly, that the wave functions 
in fact scale in a semiclassical way, ln\I' = JF/A, where depends on appropriate 
combinations of other parameters and A. This empirical evidence suggests that tun- 
neling in quantum cosmology may still be described semiclassically, but the relevant 
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classical solutions are entirely different from the instanton existing in this model. 
This conjecture is further supported by our study of the dependence on the a-0 cou- 
pling parameter (m in eq. 0) at small but fixed A. As we already mentioned, at 
small m the matter degree of freedom is not excited due to tunneling, while at larger 
m it is, and its energy behind the barrier is close to Vmax- We find numerically that 
the transition between the two regimes is very sharp in m and becomes close to step 
function as A gets smaller. This suggests that the wave function is a combination 
of two semiclassical exponentials, 

^ = Cie-^' + Cae"^^ (8) 

where the first exponential is due to tunneling along (j) = and coincides with eq. Q, 
while the second exponential comes from an entirely different classical path. In spite 
of considerable effort, we were unable to identify the classical solution responsible 
for tunneling with excitation of the matter degree of freedom, which would be of 
considerable interest to quantum cosmology. 

Overall, we have found that tunneling in quantum cosmology is quite different 
from that in quantum mechanics. This observation may have interesting implications 
to the discussion of initial conditions for the classical stage of the cosmological 
evolution. 

Appendix A contains details of our numerical procedure. We describe various 
checks of this procedure in Appendix B. These checks ensure that our results are 
reliable, even though we are dealing with very small numbers typical for tunneling. 
The calculations have been performed at the supercomputer SGI Origin 2000 of the 
Center for Computational Science of Boston University. 

2 The model and perturbative solution. 

It is convenient to rescale the variable a by introducing 

b = y/Xa 

and write the Wheeler-De Witt equation in the following form. 



>P(6) = , (9) 



where the coordinate representation in b is chosen, while the representation of op- 
erators and is not yet fixed. Here is the Hamiltonian of an oscillator with 
6-dependent frequency, 

n^ = lK + l^"(b)^" (10) 

where 

n^{b) = 00^ + M^f{b) (11) 
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and we introduced a parameter u for generalitjU. The parameters e and M are 
related to original parameters as follows, 

.2 



Ae 



m 



M' = — . 
A 



(12) 



It is worth noting that for arbitrary U{h) and eq. Q corresponds to the action 



S 



_ _ AT (U{b) - e) 



+ 



drj 



2N 2 ^ 



(13) 



where is the lapse function. 

In the original Wheeler-De Witt equation, one has 

U{h) = Ar(6/VA) = (feV2 - h^) 

and f{h) =6^. In our numerical study, however, we found it convenient to choose 
U{h) and f{h) is such a way that they are well localized. 



2.1 Modified models. 

We note that, since we are imposing free boundary conditions at 6 = (as U (0) = 0, 
/(O) = 0), the range of h can be extended to (— oo; +oo). It is then convenient to 
shift h in such a way that the maximum of the potential U{h) occurs at 6 = 0. To be 
able to impose the boundary conditions in a finite interval of 6, we will work with 
the following potential and interaction functions, 

U{h) = e-^V2 (14) 

and ^ 

gQ!i(fe-fei) _j_ ga2(6-b2) _|_ ga3(6-fe3) ■ 

We have chosen three sets of parameters ai, 02, «3 and 61, 62, &3, summarized in 
Table 1. Since not all have the same sign, the functions /i, /2 and /s all decay 
exponentially as \b\ — > 00. We have chosen values of c such that the maximum values 
of these functions are close to 1. Then the strength of h-(j) interaction is determined 
by the parameter M in eq. pi|l . 

The potential U (b) and the three interaction functions are shown in Fig. |21 The 
shapes of the interaction functions are quite different: while /i is substantial in front 
of the barrier (and hence the matter degree of reedom (p may be excited there), /2 and 
/a vanish in front of the barrier, with being non-zero in the classically forbidden 
region only. Thus, these three interaction functions constitute a representative set. 
We will see that the main features of tunneling do not depend on the shape of the 
interaction function. 

^For inhomogeneous modes, one would have = 1 + fc^, see eq. (01. Note that k is the 
conformal momentum of the scalar mode, which is related to the physical momentum as Pphys = 
k/a — k\/K/h. Even for modes with sub-Planckian physical momentum, the value of fc, and hence 
w, may be much larger than 1. 
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Table 1: Parameters of the interaction functions. 




Figure 2: The potential U{h) and the three interaction functions. 



For fixed choice of there are four parameters in the model, namely, A, e, uj 
and M. Naively, one would expect either from eq. (jHl) or eq. (|T3|l that as A — > with 
e, UJ and M fixed, tunneling through the barrier U{h) is described semiclassically, 
and, furthermore, that the matter degree of freedom makes a small effect (provided 
that tunneling occurs from a low- lying state of the 0-oscillator). Let us consider the 
corresponding perturbation theory to see whether this is indeed the case. 



2.2 Pert ur bat ive treatment. 

Let us consider the case A ^ 1 with e, u and M fixed, and treat li.^ in eq. Q 
as perturbation. We take the 0-oscillator at left of the barrier in its ground state. 
Neglecting Ti^ altogther, we obtain that in the classically forbidden region, except 
near the turning points, the dominant part of the wave function is 

where 
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and bi is the left turning point. To include the effect of perturbatively, we write 

Q-Soib) 

and obtain, to the first non-trivial order, that ^1/ obeys 

= +n4T)^^ , (16) 

where we introduced a variable r, instead of b, such that 

db 



and r = corresponds to the left turning point. Equation ()16|) is a "wrong sign" 
Euclidean Schrodinger equation; the "wrong sign" (as compared to conventional 
quantum mechanics) originates from the "wrong sign" of the Hamiltonian for the 
scale factor |2j. 

To deal with eq. (fT^. let us decompose '^^{j) into eigenstates of the "instanta- 
neous" Hamiltonian (fTUj). i.e., the eigenstates \l/^^^(r) of an oscillator of frequency 
fi(r) = n{b{r)), 

oo 

'^<p{r) = Y.Cn{r)¥^\r). 

n=0 

We find that (fT^ translates into a set of equaions for the coefficients c„, 
dcr, „ , , . 1 dVt r, T, r 1 d^ 



= Q(^n + l/2)c„ - — — ^{n + l)in + 2)cr,+2 + — — ^'^(n - l)c„_2 . (17) 

This set is still hard to solve. To get an idea of its solutions, we make use of the 
adiabatic approximation, and treat the last two terms in eq. flT7|) perturbatively. 
Since we consider an oscillator in its ground state on the left of the barrier, the 
zeroth order solution is 

T 

i fdr'nir') 

(r)=e 6^^,. (18) 

To the first order we obtain that the only non-vanishing coefficient is 



4^^(r)=e dr'^ , (19) 

The first factor in the integrand, 

1 dQ{T') _ 1 dVL(b') 



2V2^](r') dr' 2Vt{b') db' 



U{b') 



accounts for the creation of a pair of 0-quanta at point b' in the forbidden region; it 
is small in the adiabatic regime, 

1 dVLir') 



fi(r') dr' 
9 



which occurs at large uj or small M. However, there is a growing exponential factor 
in the integrand of eq. (fT^. which reflects the property that the creation of 0-quanta 
enhances tunneling. The perturbation theory about the solution ()15|1 breaks down 
when this factor becomes dominant; in that case one expects that matter has strong 
back reaction on the process of tunneling. 

Our perturbative treatment shows that at given lo and small M, tunneling in 
quantum cosmology is indeed described by the wave function (fT3j) and effects of the 
matter degree of freedom are small. However, as M increases, the tunneling regime 
changes, and matter effects become important. A similar phenomenon occurs if M 
is kept fixed and small, while lo increases instead: high energy quanta are copiously 
created even if their interaction to gravity is wealcl. This is in sharp contrast to 
conventional quantum mechanics. 



3 Numerical analysis. 
3.1 Classical solutions. 

The system with the action ()13p may either tunnel through the potential barrier or 
overcome it classically, depending on the parameters. Before studying tunneling, it is 
instructive to discuss briefly the classical solutions and find the range of parameters 
in which classical transitions over the barrier do not occur. 

Let the classical system begin its evolution from large negative h towards the 
barrier. In terms of the variables h and 

V? = VA0 

the classical equations of motion read 

d% dU{b) 1 ,df{b) 



dt^ 



n'{b)ip . (20) 



V 

Here we fixed the gauge by introducing the standard cosmic time t = J dr]'N(r]'). 
The classical constraint 

na + H^ + i = (21) 

should also be satisfied. 

Note that the parameter A does not enter the classical equations; as we pointed 
out above, it is analogous to h and appears in the quantum problem only. 

At large negative b the variables decouple. One can freely choose the value of 
b in this region at initial time. Then there are two independent initial data, which 
we choose as the 0-oscillator energy Ey^{t = 0) = E^^ and its initial phase 6^0. The 
"energy" associated with the scale factor, which we define as 

Eb = —'Ha , 



•^According to eq. 112|) this may happen even if the physical masses and momenta are small. 
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is then determined at the initial moment of time from the constraint (PT|) . 



Ebfl — E^fi + — . 
For given initial h the latter relation determines h{t = 0) 



(22) 



T^bit = 0), SO the initial 



data are completely defined. Conversely, one may choose the initial value of Tib, and 
hence Ebfl, arbitrarily, and find e from eq. ()22|) . Note that 



since e is non- negative. 

For every pair of parameters (M, u) there exists a region in the plane (-^1^,0; Ebfl) 
in which at some values of the phase ^0 classical over-barrier transitions occur. The 
rest of this plane may be called "the no-transition region" of initial data: irrespec- 
tively of the oscillator phase, the barrier cannot be overcome classically for initial 
data belonging to this region. Let us note that at E^p^ = the (/(-oscillator remains 
in its classical ground state (p = during the entire evolution, so the initial data 
{E^fi = 0; Ebfl > Umax/^ = V^) "se to classically allowed over-barrier tran- 
sitions, while the interval (-E<^,o = 0; < Ebfl < 1/A) belongs to the no-transition 
region of initial data. 

For E^pfi 7^ classical over-barrier transitions may occur even if Ebfl is lower 
than the barrier height 1/A. Indeed the parametric amplification of (/^-oscillations 
may take place and give rise to over-barrier transitions for Ebfl < 1/A. Hence, 
the boundary of the no-transition region of initial data is not a straight line; this 
boundary depends on parameters M and uj. 

We have found the regions in the plane of initial data (-E(^,o; -^6,0) in which classi- 
cal over-barrier transitions occur at some values of the initial oscillator phase 9q, and 
those where these transitions do not happen for any 6*0, by solving numerically the 
system ((201); (EH)- Details of the numerical procedure are given in Appendix A. We 
have found that even at relatively large interaction strengh M, classical transitions 
do not occur for E^^ -C 1/A and Ebfl < 1/A. An example is shown in Fig. IHl In 
what follows we always choose the parameters in such a way that the classically 
allowed over-barrier transitions are impossible. 

3.2 Solutions to the Wheeler— De Witt equation. 

Our main purpose is to study solutions to the Wheeler-de Witt equation (0) for 
our model systems numerically, without making use of any approximation. The 
tunneling boundary conditions are as follows. At large |6| the variables h and (f) 
decouple, and the frequency of the (/(-oscillator is equal to uj. Hence, the asymptotic 
wave function may be decomposed into the oscillator eigenstates. 



Ebfl > E^^Q 



(23) 



00 



(24) 



n=0 
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0.4 0.5 0.6 

A-Ei^, 

Figure 3: Boundaries of no-transition regions in the plane of initial data for u = 0.6 
and different M. The interaction function is /2. Classical over-barrier transitions do 
not occur for initial data below these lines and any initial phase of the 0-oscillator. 
The dotted line corresponds to Eb^o = E^^, see eq. (j2 



where 
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e f 1 

+ UJ \ n + 



LA V 2. 

We consider tunneling from a state of fixed oscillator excitation number uq. Then 
the tunneling boundary conditions are that the negative momentum component of 
the wave function vanishes on the right of the barrier, while the positive momentum 
component on the left of the barrier has n = no, 

= 0, 

tn = <no- (25) 

These conditions are direct analogs of the tunneling boundary conditions in quantum 
mechanics. 

Unless explicitly stated, we will consider tunneling from the ground state of the 
(/)-oscillator, 

no = . 

We will discuss the dependence on uq towards the end of this section. 

We describe our numerical procedure in Appendix A. Various checks of this 
procedure are summarized in Appendix B. Here we present the results of our calcu- 
lations. 

To describe our solutions at finite b, including classically forbidden region, it 
is convenient to work in the basis of the "instantaneous" Hamiltonian of the 0- 
oscillator. The eigenf unctions obey 
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where 

The wave function is expanded 

oo 

^{b) = Y.Cn{b)¥^\b). (26) 

71=0 

We will be interested in the behavior of the occupation numbers |C„(6) p as functions 
of the scale factor b and excitation number n. We will see that at small A, which is 
of primary interest to us, the occupation numbers are exponential, 

|a(6)poce^^(^'^^"), 

so they are sharply peaked, at given b, near a certain value of n. For sufficiently 
large values of the coupling parameter M and small A, typical plots of Aln are 
shown in Figs. ID and El 
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Figure 4: The occupation number as function of the energy of the 0-oscillator 
at different values of the scale factor, for the interaction function fi and uq = 
(tunneling from the oscillator ground state), A = 0.0225, M = 0.2, uj = 0.6, e = 0.27. 

To discuss these plots, we recall first that the potential ()14|1 is centered at 6 = 
and its width equals a/2, so the potential is small at 6 = ±2, the first and last 
points displayed in Fig. HI In Fig. El we show also the occupation numbers at larger 
b where the interaction function is large and additional creation of 0-quanta occurs, 
now in the classicaly allowed region. We also note that we consider quite small 
values of A, so the occupation numbers are actually very small (because of tunnel- 
ing suppression) in the forbidden region and behind the barrier. Nevertheless we 
are confident that the results are reliable, partially due to checks of our numerical 
procedure summarized in Appendix B. 
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It is clear from Figs. |3] and El that on the left of the barrier, the 0-oscillator is 
mostly in its ground state (since we have chosen uq = 0), while in the classically 
forbidden region its wave function gradually becomes peaked at the value of n such 
that the total matter energy, (-E^ + e/A), is close to the height of the barrier, 
Vmax = Comparing Figs. El and we see that this property is essentially 

independent of the shape of the interaction function, as well as on parameters of the 
model (provided that M is sufficiently large, see below). We have checked that this 
is indeed the case by studying the solutions for numerous sets of parameters and 
also for the interaction function /a. 

The mechanism due to which tunneling is accompanied by the creation of 0- 
quanta is clear. The (p-b interaction results in weak creation of these quanta (either 
in front of the barrier, as is the case for the interaction function /i, or in the forbidden 
region, as is the case for /2 and /s). The low modes get strongly suppressed in the 
forbidden region, while the modes of energy ~ 1/A survive. The latter effect 
occurs because the total energy of matter is close to the barrier height. Hence, there 
is strong back reaction of created 0-quanta on tunneling. 

One remark is in order. Even though we impose the boundary condition 
with no = 0, the wave function on the left of the barrier still contains a small 
admixture of the excited states of the 0-oscillator. This is due to the presence of 
the reflected wave. 

Another way of studying the properties of our solutions is to consider the average 
energy of the 0-oscillator. For finite b we define it as follows, 

g E^;\b)mb)\' 

imb) = ^ (27) 

E •ia(&)P 

n=0 
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In the asymptotic region b +00, a more appropriate definition would involve the 
flux factor, as in conventional quantum mechanics, 

(28) 

b— >+oo 

Since \Cn\ are sharply peaked at a certain n, the two deflnitions do not differ much: 
{E^y^^y^P ^ {E^){b — > +00). We will use the deflnition (j2H) for general 6, and the 
deflnition (j^Hj) when discussing the asymptotics 6 ^ +00. 



E knUj{n + l/2)\C„ 



Til kn\Cn\'^ 
n.=n 




0.2 0.4 0.6 0.8 1 



Figure 6: </)-oscillator energy at large h as function of e for the interaction function 
/i and A = 0.0256, M = 0.2, lu = 0.6. 



As one anticipates from Figs.EJandEl the total matter energy at 6 ^ +00 should 
be close to Vmax = 1/^5 

{E^yy^P ^ j{l - i) . (29) 

This is shown in Fig. IHl In fact, the relation (j29|l should not be precise, because 
the relevant quantity is the oscillator energy near the top of the barrier rather than 
its asymptotic value. For given n the former is equal to (n + l/2)Q{b ^ 0) = 
(n + l/2)^uj'^ + M'^fih ^ 0) which is larger than the asymptotic value (n + 1/2)0;. 
Hence, it is more appropriate to study the oscillator energy (j2Zj) as function of b for 
different values of parameters. The corresponding plots are shown in Fig. which 
clearly displays how the system climbs up the potential barrier, and then stays at 
constant |Cnp, with energy decreasing due to the effect just described. 

Figure [7| reveals another property: at smaller interaction strength M, the (p- 
oscillator gets excited "later", i.e., at larger b. This may be anticipated from eq. (fT^ . 
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Figure 7: Total matter energy in units of the barrier height, as function of b for 
different M; the interaction function is /s and A = 0.0225, u = 0.6, e = 0.27. 



as the exponential factor compensates for small M at larger b. In fact, the properties 
we discuss should be absent for very small M, and the perturbative regime of section 
2.2 should set in instead. We will come back to this point later on. 

It is instructive to see what happens at different A. We still consider tunneling 
from the oscillator ground state, uq = 0. At large A the barrier is effectively absent, 
so the 0-oscillator remains in its ground state with energy 



E 



(0) 



(30) 



At small A the energy of the 0-oscillator increases substantially, to 



asymp 



const 
A 



(31) 



The overall behavior of the asymptotic value of the oscillator energy is shown as 
function of A in Fig. IHlin log-log scale. One indeed observes the asymptotics (p?T|l at 
small A and ()30|1 at large A. The oscillatory behavior of {E^)"'^^"^^ as function of A, 
shown in detail in Fig. 13 may be understood as follows. For finite A, the number of 
"under-barrier" oscillator levels (which have energy below (1 — e)/A) is finite, and as 
A decreases, some over-barrier levels become under-barrier. Since the system tends 
to climb on the top of the barrier, but not much above the top, the structure of 
levels near the top matters. This structure changes with A in an oscillatory way, 
hence the oscillatory behavior shown in Fig. 121 

A similar dependence on A is inherent in the magnitude of the wave function 
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Figure 8: 0-oscillator energy at 6 ^ +oo as function of A, in log-log scale. Straight 
lines are asymptotics (jHUj) and pip. The interaction function is /a, the parameters 



are: M = 0.4, u = 0.6, e = 0.27. 
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Figure 9: The same as in Fig. |H1 but for a region of small A. The scales are linear, 
behind the barrier. Let us define a quantity F as followfl 

(32) 



r 

eA 



n=0 



b—* + OD 



where we again consider tunneling from the ground state of the 0-oscillator. We 

■^In conventional quantum mechanics, exp(r/A) would be tunneling probability, while in quan- 
tum cosmology such an interpretation may be debatable. 
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find numerically that at small A the exponent scales as 



T{M,u) 



(33) 



i.e., it is independent of A. This is shown in Fig. ^| Note that in the regime we 
consider, F/A is much larger than the naive tunneling exponent 



r da j2iV{a) - e) 

J ai 



(34) 



entering eq. (jlj). Thus, creation of 0-quanta has dramatic effect on the magnitude 
of the wave function behind the barrier. 
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Figure 10: The magnitude of the wave function behind the barrier as function of 
a/a. The quantity F is defined in eq. ((221) • Short dashed line in the lower left part 
corresponds to the naive exponent ()34|1 . The interaction function and parameters 
are the same as in Fig. [3 



In the regime of relatively large interaction parameter M, the dominant part 
of the wave function behind the barrier corresponds to (^-oscillator excited to a 
particular energy, such that the total matter energy near the top of the barrier is 
close to Vjnax = 1/A. Heucc, this state is essentially independent of the state of the 
0-oscillator on the left of the barrier, i.e., of no, up to the overall penetration factor 
(the latter depends on uq strongly). This is illustrated in Fig. HU 

Let us now discuss the case of small M, when the couping between and b 
is weak. According to the perturbative analysis of section 2.2 in that case the 
0-oscillator does not get excited due to tunneling, provided that u is small. How- 
ever, one may expect that at large enough u the situation is different and that 
the system is back in the regime of strong particle creation. Hence, as uj in- 
creases, the asymptotic energy of the oscillator should change from a small value to 
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Figure 11: 0-oscillator energy at 6 — +oo as function of the oscillator energy on 
the left of the barrier for interaction function /a and A = 0.0324, M = 0.2, e = 0.27. 
The linear behavior at large Ea^q corresponds to classically allowed transitions. 



l^E^Y'^y^'P K, (1 — e)/A. This is shown in Fig. ^| Of course, the value of uj at which 
the transition between the two regimes occurs depends on M and other parameters, 
but in any case the scalar quanta with large enough oj are copiously created due to 
tunneling. In the context of quantum cosmology this means that the Universe after 
tunneling is filled with particles, possibly of large mass and high spatial momenta. 



3.3 Scaling with A. 

We have seen in section 3.2 that at small A and sufficiently large M and/or cu, the 
</)-oscillator gets excited in such a way that the total energy of matter is close to the 
height of the barrier. This means that the wave function of the oscillator is peaked 
near a certain excitation number n oc 1/A. One may consider this property as an 
indication that the tunneling transition is of semiclassical type. Further support of 
this conjecture comes from the semiclassical-type scaling of the overall magnitude 
of the wave function behind the barrier, eq. ()33j) . Let us describe further empirical 
evidence in favor of the semiclassical nature of the tunneling transitions in the regime 
where the 0-oscillator gets strongly excited. 

One observation is that not only the overall magnitude, but also the entire wave 
function has the semiclassical form 

C^(h) = e*^(^-^<') , (35) 

both in the classically forbidden region and behind the barrier. As an example, we 
show in Figs. El and El real and imaginary parts of JF = A lnC„ behind the barrier 
as functions of AE^^ at different values of A. Scaling property ()35|) is clearly seen in 
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Figure 12: 0-oscillator energy at 6 ^ +oo as function of uj for the interaction 
function /2, M = 0.3, A = 0.04, e = 0.27 and Uq = 0. At small u the oscillator does 
not get excited, while at larger u it is strongly excited due to tunneling. Oscillations 
at large u have the same origin as oscillations in Figs. |Hl El and ^1 



these figures. We observed very similar scaling of the wave function in the classically 
forbidden region, though its shape is quite different there (the latter is shown for a 
representative value of A in Figs. |3] and Ej). 
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Figure 13: Real part of the exponent of the oscillator wave function behind the 
barrier at different A for interaction function fi and M = 0.2, u = 0.6, e = 0.27. 
Note that the wave functions themselves differ by more than 10 orders of magnitude 
at different values of A. 
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More evidence comes from the study of the transition from the perturbative 
regime, occuring at small values of the coupling parameter M, to the regime of 
strong excitation of the 0-oscillator at larger M. As shown in Fig. this transition 
becomes more sharp as A decreases, and becomes essentially step function at small 
A. We consider this as an indication that the wave function behind the barrier and 
in the forbidden region is a combination of two exponentials of semiclassical type, 
eq. ©. 

What classical solutions describe the tunneling transitions with strong excitation 
of the 0-oscillator, remains an open problem. 
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Figure 15: 0-oscillator energy behind the barrier as function of M at different A 
for the interaction function fi and u = 0.6, e = 0.27. 

A.l Classical solutions. 

To obtain the solutions to eqs. (pUjl and (jHIJ, we chose the initial value b = —15, well 
outside the barrier, and imposed the constraint (PT|) at the initial moment of time. 
We found the classical solution by making use of the fourth order Runge-Kutta 
method with adaptive step size (see, e.g., Refs. [12]) for the first order Hamiltonian 
equations of motion. We checked that the constraint (j25) was satisfied during entire 
evolution with relative precision 10~^. For given initial oscillator phase ^o, the 
boundary of the no-transition region of initial data was found with precision 10~^ 
by "division by 2" search in Eh^. The phase Oq was discretized into 100 points, and 
the minimum in Oq of the boundary of the no-transition regions was found. 

A. 2 The discretized Wheeler— De Witt equation. 

We solved the Wheeler-De Witt equation (0) by expanding the wave function in 
eigenstates of the oscillator of fixed frequency iv, 



where "^^^ are independent of bu Asymptotically, Cn coincide with the coefficients 
entering eq. (j21|), so the boundary conditions are precisely ((23) • The Wheeler- 
De Witt equation takes the form of an ordinary differential equation for the vector 



^ This expansion, which is more convenient for computational purposes, is related to the ex- 
pansion into "instantaneous" modes of eq. H26|l . see A. 5. 



n 




(36) 
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where the matrix A is defined as 



1 / M'^f(b)' 

-(f/(6)-e)-^(n + l/2) 1 + 



^n' n 



+ 



2uj 



{j{n + l){n + 2)5n', n+2 + \ln{n - „_2^ . 



To solve eq. numerically, we truncated the system to a finite number of levels, 
i.e., set < n < (A^^o — 1). We considered the system in a finite interval of h and 
discretized the variable h by introducing a lattice of spacing A. The sites were chosen 
at points 

6j = iA , 

where i = —Nb, . . . Nb. 

One way to discretize eq. (j36|) would be to make use of the symmetric approxi- 
mation for the second derivative 



d^Cjb) 



1 

A2 



[C{b + A) + C{b - A) - 2C{b)] + 0{A' 



(37) 



However, we made use of the Numerov-Cowling method which has higher precision. 
Namely, the Taylor expansion of C{b) gives 



C'i+i + C'i-i ~ 2(7^ — A^ 



dh^ 



+ 



A^ 
12 



d^C 



+ 0(A« 



Making use of eq. (jHUj) . one writes the second term on the right hand side as 



UAdl^ 



(AC) 



A^ 
12A 



(3^ 



where eq. fl37|) was used for the derivative of the vector (AC). In this way we obtain 
the discretized version of eq. ()36|) with error O(A^), 



A2 \ ~ A' 
12A / I 12A 



-A 



a 



0. 



(39) 



The unknowns of this finite difference equation are numbers Cj^„ with i = —Nb, . . . Nb, 
n = 0, . . . {Nq — 1). So, the number of unknowns is No{2Nb + 1). The system (jH^jl 
is a set of N(){2Nb — 1) equations. The boundary conditions (j2Sl) after discretization 
become 2No linear equations 



^n, no ^1 



2ik„A 



Nb, n 







(40) 



Hence, the total number of equations coincides with the number of unknowns, so 
Ci^n are uniquely determined. 

As we will see below, the solution of the discretized system well approximates 
the solution of the original one only if the numbers of lattice points and oscillator 
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levels are quite large, A^^ ~ 20000, Nq ~ 150. Numerical methods for generic system 
of equations would fail in this case, since the total number of complex equations 
is more than 3 million. However, following Ref. we exploited the special form 
of eq. to derive the solution by an iterative procedure. Indeed, by inverting 
numerically a set of A^o x matrice^ one writes eq. in the following form, 

Ci = LiCi-i + RiCi^i , (41) 

where Li and Ri are again A^o x matrices. Equations (j^T|) can now be used to elim- 
inate a set of Ci variables with given i, leading to new equations, which, with some 
matrix algebra, can be cast again into the form (PT|) . Thus one can progressively 
eliminate all the variables Ci at the intermediate points i = {—Ni, + 1) . . . {Ni, — 1) 
expressing them in terms of linear combinations of C-n^ and Cn^- By substituting 
these linear combinations into eq. ()40j) . we obtained a system of linear inhomoge- 
neous equations for the complex vectors C^n^ and C^^ which was solved in straight- 
forward way. Then the actual values of C-n^ and Cjy^ were introduced back into 
expressions for the rest of unknowns, and in this way the complete wave function 
was found. 

It is worth pointing out that the procedure in which intermediate variables are 
excluded involves operations with real matrices only. Note also that the interaction 
between different C*„ in eq. (jHUj) occurs only if they have the same parity. Hence, it 
is sufficient to consider only even (or only odd) modes. 

The algorithm described above lends itself to parallel implementation: the vari- 
ables Ci at odd sites are iliminated in parallel, and the procedure is repeated in a 
recursive manner. Our code scaled well with the number of processors A^proc, and 
most of the results presented in this paper were obtained with A^proc = 64. 

One must of course be careful that the effects of discretization and truncation 
of the original Wheeler- De Witt equation are kept small. We used a value for 
the cut-off of the number of modes A"o such that the truncation energy exceeds the 
height of the barrier by a factor of 2, 



The lattice spacing was chosen equal to 1/12 of the minumum wavelength along b, 

1 27r V^irA 

A 



12 kNo-i 12V2Te 
We used the lattice of the size 

so that the potential is very small near the end points. We performed various 
checks, summarized in Appendix B, which show that this choice of parameters is 
appropriate. 



^An additional simplification is that the matrices (2 — 5Ai/6) that must be inverted are sparse, 
with non-zero entries on three diagonals only. Inverting matrices of this type takes much less 
computer time than with to matrices of generic form. 
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A. 3 Iterative improvement of precision. 

In our calculation we used a double precision representation of real numbers, which 
entails round-off errors of order 10~^^ to 10~^^. The exponential behavior in b and n 
gives origin to numbers Cj^„ varying over a range that extends well beyond the limits 
of double precision. What makes it still possible for us to use double precision is 
that only coefficients C*j^„ with neighboring values of i and n and comparable magni- 
tude are coupled by the equations that must be solved. Thus one never encounters 
additions or subtractions of variables with wildly differing exponents. The elimina- 
tion procedure used for solving the equations can nevertheless propagate round-off 
errors, so that special care is needed to ensure that the values of the smaller Cj_„ are 
accurate. We accomplished this by an iterative improvement of precision, which we 
briefly describe here. 

Let C^^n be the initial solution obtained by the elimination procedure. The 
propagation of round-off errors will give origin to absolute errors of order e = 10~^^ 
to 10~^^ in the coefficients Cj^„. Thus those coefficients whose true value is smaller 
than e will turn out to be totally wrong and the equations coupling them will have 
errors of the same order of magnitude. It is convenient to denote these equations 
symbolically by AC = 0. Substituting into the equations the numbers C^^^ we 
obtain instead AC^^^ = B^'^\ where B^^^ ~ e. We introduce therefore a correction 

which we find by solving A SC^°^ = -B^°\ Note that the equations for 
only involve right-hand side terms of order e, so that the propagation of round- 
off errors will now give origin to errors of much smaller absolute magnitude ~ e^. 
With (5C(o) we form the first iterate of the solution C(i) = C^") + 6C'^°\ While the 
numbers C^^^ and 5C^^^ will typically not be smaller than e, the first iterates C^^^ 
corresponding to coefficients whose true value is less than e will accordingly take 
smaller values, being now affected by errors ~ e^. Inserting the first iterate C-^^ 
into the equations we find now errors B^^'' = A C^^^ ~ and the procedure can be 
repeated until the values of the smallest coefficients become stable. 

We show in Fig. El the relative error as function of n after different numbers of 
iterations. This quantity is defined as the ratio of the modulus of the left hand side 
of eq. (PT|) to the sum of moduli of all terms in this equation. The absolute value 
of Ci^n at 6 = at different steps of the iteration procedure is presented in Fig. UHl 
We conclude that the iteration procedure indeed enables one to obtain the solutions 
with relative, not absolute, error of order 10~^^. 

A. 4 Total computational time. 

Our main emphasis in this paper is the study of the solutions to the Wheeler-De Witt 
equation at small A. Let us estimate the dependence of the time of calculations on 
this parameter. At each step of the iteration procedure of section A. 3, A^^ ■ A'q 
floating point operations are performed. The number of iterations is proportional to 
logarithm of the minimum value of the wave function. The latter scales as A^j^g^- oc 
1/A, see eq. Thus, the time of calculations is proportional to 

Nt ■ ■ A^iter « ^ • 
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Figure 16: Wave function at 6 = at different steps of the iteration procedure. 
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Figure 17: Relative error at different steps of the iteration procedure. 

This estimate shows that reaching lower values of A occurs with dramatic increase 
of processor time. 

In practice, for A = 0.0225 and uo = 0.6, the time of calculations with 64 proces- 
sors was equal to 30 to 40 minutes. At larger values of u the time of calculations 
was somewhat smaller. 
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A. 5 Bogoliubov transformation. 

Numerical calculations were conveniently performed in the basis of the eigenfunc- 
tions of the oscillator of 6-independent frequency uj. The interpretation of the results 
is done, however, in the basis of the "instantaneous" Hamiltonian, describing the 
oscillator of frequency So, for every h we performed the Bogolubov transfor- 

mation and obtained the coefficients C„ starting from Cn- 

The key element of this procedure is a set of transformation coefficients 

a" 1 ^^""^ 

Here ccq denotes the annihilation operator in the basis of the "instantaneous" Hamil- 
tonian. The first factor on the right hand side is given explicitly by 

while the matrix elements {ki_^\cYQ\mi_^) can be calculated numerically, using the har- 
monic oscillator algebra. 



B Checks of numerical procedure. 

Since we were dealing with very small values of the coefficients Cj,„, we had to 
perform a series of checks of the entire procedure. To make sure that round-off errors 
do not spoil our calculations, we repeated the calculations with quadruple precision 
(32 decimal places), and found that the results for the wave function coincided with 
those obtained with double precision, with relative error 10~^^. This check was 
performed for the interaction function /s and the following values of parameters: 
VA = 0.3, 0.21 and 0.17; M = 0.2; u = 0.6; e = 0.27. 

We also checked that the effects of discretization of b and truncation of the 
oscillator energy levels are small, by varying the parameters Nj,, Nq and A by a 
factor of 2. The interaction function in this check was /a, the parameters were 
= 0.21; M = 0.2; cu = 0.6; e = 0.27. The relative difference between the 
solutions was always less than 10~^. 

Another way to check our procedure is to make use of the current 

J = 1 (c+dbC - dbC+C) , 

which is conserved in the continuum theory as a consequence of the Wheeler-De Witt 
equation, i.e., 

while the discretized version does not have this property. By calculating the value 
of the current for every lattice site, we checked the discretization error. The pre- 
cision of the calculation of J for given set of vectors Ci must be O(A^), otherwise 
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the conservation of current would hold with lower precision than the precision the 
solution was found. From the Taylor expansion we find 

30 ^TT^' + 40 ^\~T^* + 72 THC^'-' + 

AV+(5.'^).^ A^+^m^ +^r+'^^r + r42) 

20 ^^^^+1 - ^^'^+1 + ^^^^ 

+ 360^^ nc^'^' + 360^^ 

Here we made use of eq. p6|) to express higher derivatives of Cj. The first and 
second derivatives of the matrix A were found by straightforward differentiation of 
the expression (p?7j) . 

Making use of the calculated current J, we estimated the relative error of the 
calculation of 

^7Vo-l|C*jp 

where the exact value of the current was found in the asyptotic region of large 
positive 6, 

J,,act = lim E '^kn\Cn{h)? ^ E 2A;„|C'^„n|l (44) 

\n=0 / n=0 

Because of the exponential suppression of the wave function behind the barrier, 
eq. pij) is satisfied with better precision, as compared to the precision at which the 
current is calculated in the classically forbidden region. Equations (HH), PHj) and 
()44|1 give an estimate of the relative error in the calculation of Cj,„, which we found 
to be less than 10~^. 

As one more check we changed the sign of the Hamiltonian of the scale factor, 
and obtained a conventional quantum mechanical system of two degrees of freedom. 
We checked for that system, that the sum of reflection and transmission coefficients 
was equal to 1 with precision 10~^°. 

Finally, we compared the numerical solutions at small M and A to the analytic 
solutions obtained within the adiabatic approximation of section 2.2. The first two 
modes of these two solutions coincided within 3 %. The latter, fairly large deviation 
we ascribe to the approximations involved in obtaining the analytic solutions. 
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